subroutine readsink(lunit,t_array0,myrank)
  use param
  use data_param,only:gsfn
  use data_param,only:geom_exist
  use data_param,only:sink_exist
  use data_param,only:n_sink_data
  use label_data_param,only:n_obj_type
  use geom_data_param,only:n_geom_obj
  use sink_data_param,only:sink_map
  use sink_data_param,only:sink_type
  use sink_data_param,only:n_sink_type
  implicit none

  interface
     function reallocSinkReactType(p,n)
       use derived
       type(SinkReactType),pointer,dimension(:) :: &
            reallocSinkReactType,p
       integer,intent(in) :: n
     end function reallocSinkReactType
     
     function reallocReal(p,n)
       use param
       real(dbl),pointer,dimension(:) :: reallocReal,p
       integer,intent(in) :: n
     end function reallocReal

     function reallocInt(p,n)
       integer,pointer,dimension(:) :: reallocInt,p
       integer,intent(in) :: n
     end function reallocInt

     integer function str2label(str)
       character :: str*(*)
     end function str2label

     function reallocSinkProductType(p,n,q)
       use derived
       type(ProductType),pointer,dimension(:) :: &
            reallocSinkProductType,p
       integer,intent(in) :: n,q
     end function reallocSinkProductType
  end interface

  integer,intent(in) :: lunit,t_array0(8),myrank

  integer :: nunit,count_line,len_str,asci_code
  character(len=len_linebuf) :: eatline
  character(len=len_wordbuf),dimension(max_n_word) :: eatword
  character(len=8) :: d_time_str

  integer :: i
  integer :: gid,lid,sid,ns,nd,np

  if (.not.sink_exist) return
  if (.not.geom_exist) return
  if (n_geom_obj.eq.0) return
  
  nunit=1

  open(unit=nunit,file=gsfn,status='old',err=103)

  count_line=0
  ns=0
  do
     read(nunit,'(a)',end=100,err=101) eatline
     count_line=count_line+1
     len_str=len_trim(eatline)
     if (len_str.eq.0) cycle
     eatline=adjustl(eatline)
     asci_code=ichar(eatline(1:1))
     if (asci_code.eq.ccode) cycle
     read(eatline(1:len_str),*,end=102) &
          (eatword(i),i=1,max_n_word)
102  continue

     ! c1:geom_obj
     ! c2:obj_type
     ! c3:energy barrrier
     ! c4:prefactor
     ! c5:interaction distance
     ! c6:number of products
     ! c7-:types of products

     ! read the index of geometric object
     read(eatword(1),*) gid
     if (gid.lt.1.or.gid.gt.n_geom_obj) then
        call terminal_msg1(trim(adjustl(eatline)))
        call abnormal_stop("readsink",&
             "The index of geometric object is out of range.")
     end if

     ! read the label of the reactant
     lid=str2label(eatword(2))
     if (lid.eq.0) then
        call terminal_msg1(trim(adjustl(eatline)))
        call abnormal_stop("readsink",&
             "Unrecognized label is found in sink file.")
     end if

     ! Count the number of data
     n_sink_data=n_sink_data+1

     ! create geometric reaction type
     sid=sink_map(gid,lid)
     if (sid.eq.0) then
        ! create sink type
        ns=ns+1
        n_sink_type=ns
        sink_type=>reallocSinkReactType(sink_type,ns)
        sink_map(gid,lid)=ns
        sink_type(ns)%gid=gid
        sink_type(ns)%lid=lid
        sink_type(ns)%n_data=1
        sid=ns
     else
        ! add to the existing list
        sink_type(sid)%n_data=&
             sink_type(sid)%n_data+1
     end if

     ! Number of data in this sink type
     nd=sink_type(sid)%n_data

     ! Re-allocate the arrays for reaction data
     sink_type(sid)%react_eb=>reallocReal(&
          sink_type(sid)%react_eb,nd)
     sink_type(sid)%react_nu=>reallocReal(&
          sink_type(sid)%react_nu,nd)
     sink_type(sid)%react_rd=>reallocReal(&
          sink_type(sid)%react_rd,nd)
     sink_type(sid)%react_rd2=>reallocReal(&
          sink_type(sid)%react_rd2,nd)
     sink_type(sid)%react_const_rate=>reallocReal(&
          sink_type(sid)%react_const_rate,nd)
     sink_type(sid)%n_new=>reallocInt(&
          sink_type(sid)%n_new,nd)
     sink_type(sid)%new=>reallocSinkProductType(&
          sink_type(sid)%new,nd,sid)

     ! read the data
     read(eatword(3),*) sink_type(sid)%react_eb(nd)
     read(eatword(4),*) sink_type(sid)%react_nu(nd)
     read(eatword(5),*) sink_type(sid)%react_rd(nd)
     sink_type(sid)%react_rd2(nd)=&
          sink_type(sid)%react_rd(nd)*&
          sink_type(sid)%react_rd(nd)
     read(eatword(6),*) sink_type(sid)%n_new(nd)

     np=sink_type(sid)%n_new(nd)

     ! annihilation
     if (np.eq.0) cycle

     ! types of reactants
     allocate(sink_type(sid)%new(nd)%prod_type(np))

     do i=1,np
        lid=str2label(eatword(6+i))
        if (lid.eq.0) then
           call terminal_msg1(trim(adjustl(eatline)))
           call abnormal_stop("readsink",&
                "Unrecognized label is found in sink data.")
        end if
        sink_type(sid)%new(nd)%prod_type(i)=lid
     end do

  end do

103 continue
  call abnormal_stop("readgeom",&
       "Error occurs while the geometric file opens.")

101 continue
  call abnormal_stop("readgeom",&
       "Error occurs during reading the geometric file.")
  
100 continue

  close(nunit)

  ! Print the message to the log file
  call log_msg1(lunit,"")
  call log_msg1(lunit,"Message from the subourinte (readsink)")
  call log_msg1(lunit,"")
  write(lunit,'(1x,i0,1x,"lines read from the ",a)') &
       count_line,trim(gsfn)
  write(lunit,'(1x,"Number of types for sinks=",1x,i0)') &
       n_sink_type

  call print_d_time(t_array0,d_time_str)
  if (myrank.eq.0) &
       write(*,'(a8,":READSINK DONE.")') d_time_str(1:8)
  write(lunit,*)
  write(lunit,'(a8,":READSINK DONE.")') d_time_str(1:8)  

  return

end subroutine readsink
